#include<set>
#include<vector>
#include<map>
#include<math.h>
#include<iostream>
#include<fstream>
#include<sstream>
#include<cstdlib>
#include<gsl/gsl_randist.h>

#ifdef USE_GPU
#include<CL/cl.hpp>
#include"../common/clsafe.h"
#endif
#include"hmm_cnv_dimension.h"
#include"../common/main.hpp"
#include"../common/analyzer.hpp"
#include"../common/io.hpp"
#include"../common/utility.hpp"
#include "hmm_cnv.hpp"
#include "compute_alpha.hpp"

using namespace std;


void HMM_CNV::init(const ptree & pt){
  string t =  pt.get<string>("use_opencl");
  use_gpu = !t.compare("true")?true:false;
  platform_id = pt.get<int>("opencl.platform_id");
  platform_type = pt.get<string>("opencl.platform_type");
  device_id = pt.get<int>("opencl.device_id");
  kernel_path = pt.get<string>("opencl.kernel_path");
  observed_filelist = pt.get<string>("observed_filelist");
  total_markers = pt.get<int>("total_markers");
  sd_lrr_het = new float[5];
  sd_lrr_loh = new float[5];
  sd_baf_percent = new float[100]; 
  memset(sd_lrr_het,0,sizeof(float)*5);
  memset(sd_lrr_loh,0,sizeof(float)*5);
  memset(sd_baf_percent,0,sizeof(float)*100);
  for(int i=0;i<=4;++i){
    ostringstream oss;
    oss<<"lrr_sd_cn.het."<<i;
    sd_lrr_het[i]=pt.get<float>(oss.str());
  }
  for(int i=2;i<=2;++i){
    ostringstream oss;
    oss<<"lrr_sd_cn.loh."<<i;
    sd_lrr_loh[i]=pt.get<float>(oss.str());
  }
  int percents[]={0,25,33,50};
  for(int i=0;i<4;++i){
    ostringstream oss;
    oss<<"baf_sd_percent."<<percents[i];
    sd_baf_percent[percents[i]]=pt.get<float>(oss.str());
  }
  sd_lrr_normal = pt.get<float>("lrr_sd_normal");
  sd_baf_homo_del = pt.get<float>("baf_sd_homo_del");
  sd_alpha=pt.get<float>("alpha_sd");
  normal_cdf = math->normal_cdf(0,.5,sd_baf_homo_del);

  this->math = new MathUtils(123);
  // initialize the observed data
  this->pos_vec = new int[total_markers];
  //this->alpha_block_vec = new int[total_markers];
  this->rslist = new string[total_markers];
  this->gap_vec = new int[total_markers];

  cerr<<"Loading observed files\n";
  if (!load_observed_files(observed_filelist)){
    cerr<<"Could not load the set of input files\n";
    exit(1);
  }

  //set_global_alpha();
  //global_alpha = 1;

  // This method will determine the number of states
  this->states = 15;
  cerr<<"Total States: "<<states<<endl;
  state_mu_lrr = new float[states];
  state_sd_lrr = new float[states];
  sample_mu_lrr_adjusted = new float[samples];
  for(int sample=0;sample<samples;++sample){
    sample_mu_lrr_adjusted[sample] = 0;
  }
  sample_likelihood_log = new float[samples];
  state_het_ratio = new float[states];
  state_total_alpha = new float[states];
  state_tumorcn_vec = new int[states];
  state_het_ident = new bool[states];
  state_alpha_ident = new bool[states];
  state_alpha_prior = new float[states];

  hmm_matrix_size = states * total_markers;
  trans_matrix_size = states * states;
  
  this->sample_trans_mat = new float[samples*trans_matrix_size];
  this->trans_mat = new float[trans_matrix_size];
  //cout<<"Transition matrix:\n";
  for(int s=0;s<samples;++s){
    for(int i=0;i<states;++i){
      // initialize transition matrix
      for(int j=0;j<states;++j){
        //trans_mat[j*states+i] = 1./states;
        sample_trans_mat[s*trans_matrix_size+j*states+i] = i==j?.99:HMM_SMALL/(states-1.);
      }
    }
  }

  this->emission_cache_mat = new float[samples*hmm_matrix_size];
  this->forward_mat_log = new float[samples*hmm_matrix_size];
  this->backward_mat_log = new float[samples*hmm_matrix_size];
  this->prob_mat_log = new float[samples*hmm_matrix_size];
  cerr<<"Allocated large matrices\n";
  if(use_gpu){
    init_gpu();
  }
  cerr<<"Done init\n";
/**
  float tempmat[states*states];
  scaleTransitionMat(tempmat, 123);
  for(int i=0;i<states;++i){
  float rowsum=0;
  for(int j=0;j<states;++j){
  cerr<< " "<< tempmat[i*states+j];
  rowsum+=tempmat[i*states+j];
  }
  cerr<<" "<<rowsum<<endl;
   }
**/
};

HMM_CNV::HMM_CNV(){
}

bool HMM_CNV::load_observed_files(const string & filelistpath){
  ifstream ifslist(filelistpath.data());
  if (!ifslist.is_open()){
    cerr<<"File list is not found.\n";
    return false;
  }
  filelist_vec.clear();
  string line;
  // skip header
  getline(ifslist,line);
  float epsilon=.01;
  while(getline(ifslist,line)){
    istringstream iss(line);
    string file,annot;
    float g_alpha;
    iss>>file>>annot>>g_alpha;
    if (g_alpha==0) g_alpha+=epsilon;
    else if (g_alpha==1) g_alpha-=epsilon;
    filelist_vec.push_back(file);
    sample_annotations.push_back(annot);
    global_alphas.push_back(g_alpha);
  }
  ifslist.close();
  samples = filelist_vec.size();
  if (samples==0) {
    cerr<<"Must have at least one sample\n";
    return false;
  }
  this->lrr_vec = new float[samples*total_markers];
  this->baf_vec = new float[samples*total_markers];
  for(int i=0;i<samples;++i){
    ifstream ifsobs(filelist_vec[i].data());
    if (!ifsobs.is_open()){
      cerr<<"File "<<filelist_vec[i]<<" is not found.\n";
      return false;
    }
    int lastpos = 0;
    for(int j=0;j<total_markers;++j){
      getline(ifsobs,line);
      istringstream iss(line);
      iss>>rslist[j]>>chr>>pos_vec[j]>>baf_vec[i*total_markers+j]>>lrr_vec[i*total_markers+j];
      gap_vec[j] = pos_vec[j]-lastpos;
      lastpos = pos_vec[j];
    }
    ifsobs.close();
    cerr<<"Successfully loaded "<<filelist_vec[i].data()<<endl;
  }
  return true;
}

void HMM_CNV::run(){
  for(current_sample=0;current_sample<samples;++current_sample){
    init_state_objects(current_sample);
    cerr<<"\n\nBeginning HMM_CNV on sample "<<current_sample<<" with stromal contamination of "<<global_alpha<<endl;
    int sd_start = 1;
    int sd_end = 1;
    for(int sd_inflate=sd_start;sd_inflate<=sd_end;sd_inflate+=4){
      cerr<<"SD inflation = "<<sd_inflate<<endl;
      sd_baf_inflation=sd_inflate;
      sd_lrr_inflation=sd_inflate;
      compute_emission_cache_mat(false);
      sample_mu_lrr_adjusted[current_sample] = get_lrr_offset();
      cerr<<"Adjusted LRR mean for normal state is now: "<<state_mu_lrr[normal_state] + sample_mu_lrr_adjusted[current_sample] <<endl;
      compute_emission_cache_mat(true);
      // before baum welch training, initialize trans mat
      for(int i=0;i<states;++i){
        for(int j=0;j<states;++j){
          trans_mat[i*states+j] = sample_trans_mat[current_sample*
          trans_matrix_size +i*states+j];
        }
      }
      if(this->use_gpu){
      }else{
        this->forwardbackward_cpu();
      }
      for(iteration=0;iteration<3;++iteration){
        if(this->use_gpu){
    //      this->forwardbackward_gpu();
          //this->baumwelch_gpu(); // mle of transition matrix
        }else{
          this->baumwelch_cpu();
          this->forwardbackward_cpu();
        }
      }
      // after baum welch training, save trans mat
      for(int i=0;i<states;++i){
        for(int j=0;j<states;++j){
          sample_trans_mat[current_sample*trans_matrix_size +i*states+j] = 
          trans_mat[i*states+j];
        }
      }
      marginalize_posteriors();
    }
  }
}

//void HMM_CNV::set_global_alpha(){
//  float sd = .15;
//  float increment = (max_alpha-max_alpha)/(1.*steps_alpha);
//  map<float,int> alpha_counts;
//  for(float alpha=min_alpha;alpha<=max_alpha;alpha+=increment) alpha_counts[alpha] = 0;
//  for(unsigned int i=0;i<bins.size();++i){
//    bin_t & bin = bins.at(i);
//    float maxLL = -1e20;
//    float maxAlpha = -1;
//    for(float alpha=min_alpha;alpha<=max_alpha;alpha+=increment){
//      float LL=0;
//      for(int j=0;j<bin.size();++j){
//        feature_t feature = bin.get_feature(j);
//        float like_homo = math->normal_pdf(feature.maf,0,sd); 
//        float like_het = math->normal_pdf(feature.maf,(1-alpha)*.5,sd); 
//        LL+=log((1-HET_RATIO)*like_homo+HET_RATIO*like_het);
//      }
//      //cerr<<"Bin: "<<i<<", Alpha: "<<alpha<<", LL: "<<LL<<endl;
//      if (LL>maxLL){
//        maxLL=LL;
//        maxAlpha=alpha;
//      }
//    }
//    //cerr<<"maxAlpha: "<<maxAlpha<<endl;
//    alpha_counts[maxAlpha]+=1;
//  }
//  float maxAlpha = -1.;
//  int maxCounts = -1;
//  for(map<float,int>::iterator it = alpha_counts.begin();it!=alpha_counts.end();it++){
//    //cerr<<"ALPHA COUNT FOR ALPHA "<<it->first<<" IS "<<alpha_counts[it->first]<<endl;
//    if (alpha_counts[it->first]>maxCounts){
//      maxAlpha = it->first;
//      maxCounts = alpha_counts[it->first];
//      //cerr<<"MAX ALPHA IS NOW "<<maxAlpha<<endl;
//    }
//  }
//  global_alpha = maxAlpha;
//  cerr<<"Global alpha estimate at: "<<global_alpha<<endl;
//}

HMM_CNV::~HMM_CNV(){
	cerr<<"Destructor\n";
	//delete[] forward_mat;
	//delete[] backward_mat;
	//delete[] prob_mat;
}


void HMM_CNV::scaleTransitionMat(float * scaled_mat,int gapint){
        float gap = 1.*gapint;
	for(int i=0;i<states*states;++i) scaled_mat[i] = trans_mat[i];
        return;
        bool debug = false;
        
        ofstream ofs_scaled;
        if (debug) ofs_scaled.open("debugging/scaled.cpu");
        //ofs<<"Gap: "<<gap<<endl;
        //float offdiag = 1.-HMM_TRANSITION_DIAG;
	for(int i=0;i<states;++i){
		float rowsum=0;
		for(int j=0;j<states;++j){
                 	if (i!=j){
				scaled_mat[i*states+j] = trans_mat[i*states+j]*(1-exp(-gap/HMM_TRANSITION_LARGE));
				//scaled_mat[i*states+j] = trans_mat[i*states+j]*(penalty_mat[i*states+j])*(1-exp(-gap/HMM_TRANSITION_LARGE));
				//scaled_mat[i*states+j] = trans_mat[i*states+j]*(1-(penalty_mat[i*states+j])*exp(-gap/HMM_TRANSITION_LARGE));
				rowsum+=scaled_mat[i*states+j];
                        }else{
                           scaled_mat[i*states+j] = trans_mat[i*states+j];
                        }
		}
                //cerr<<"rowsum: "<<rowsum<<" "<<offdiag<<endl;              
                scaled_mat[i*states+i] = 1.-rowsum>small?(1.-rowsum):scaled_mat[i*states+i];
          //      scaled_mat[i*states+i] *= rowsum>offdiag?rowsum/offdiag:1.;
                //scaled_mat[i*states+i] = rowsum>trans_mat[i*states+i]?rowsum:trans_mat[i*states+i];
         //       rowsum+=scaled_mat[i*states+i];
		//cerr<<"rowsum: "<<rowsum<<" inflated "<<inflated<<endl;
		for(int j=0;j<states;++j){
	//		scaled_mat[i*states+j] /= rowsum;
                    if (debug) {
                      if (j) ofs_scaled<<" ";
                      ofs_scaled<<j<<":"<<scaled_mat[i*states+j];  
                    }
		}
                if (debug) ofs_scaled<<endl;
	}
        if (debug) {
          ofs_scaled.close();
          exit(0);
        }
}

float HMM_CNV::getBinomSum(int n, float pop_b, float baf){
	float binom_sum = 0;
	float sd = .1;
	float mean=.5;
	float mean_increment = 0;
	if (n>1){
		mean = 0;
		mean_increment =1./(n-1);
	}
	for(int k=0;k<n;++k){
    //float a = math->binom_prob(n-1,k,pop_b);
    //float b = math->normal_pdf(baf,mean,sd);
		binom_sum+=math->binom_prob(n,k,pop_b) * math->normal_pdf(baf,mean,sd);
// cerr<<"genotype: "<<k<<" Mean : "<<mean<<" binom: "<<a<<" normal: "<<b<<endl;
		mean+=mean_increment;
	}
        //cerr<<"binomsum: "<<binom_sum<<endl;
	return binom_sum;
}



float ordered_binom(int n,int k,float p){
  float prob = 1;
  for(int i=0;i<k;++i){
    prob*=p;
  }
  for(int i=0;i<n-k;++i){
    prob*=(1-p);
  }
  return prob;
}

// this needs to be called several times for each state

float HMM_CNV::get_baf_mean(int bac, int cn,float alpha){
  if(cn>=2 && bac==0 || bac==cn){
    // if tumor homo , assume normal is always homo
    return bac/cn;
  }
  int normal_bac = 1;
  int normal_cn = 2;
  float mu = ((1.-global_alpha-alpha)*bac + (global_alpha+alpha)*normal_bac)
      / ((1.-global_alpha-alpha)*cn+(global_alpha+alpha)*normal_cn);
  //cerr<<"alpha "<<alpha<<" tumor bac,cn "<<bac<<","<<cn<<" bafmu: "<<mu<<endl;
  return mu;
}


float HMM_CNV::get_baf_sd(int state, int bac){
  //float & alpha = state_alpha[state];
  float & total_alpha = state_total_alpha[state];
  // the variance of the normal tissue
  int tumorcn = state_tumorcn_vec[state];
  //cerr<<"\nFetched tumor CN "<<tumorcn<<" with het "<<state_het_ratio[state]<<" and BAC "<<bac<<endl;
  float frac = 1.*bac/tumorcn;
  //cerr<<"frac "<<frac<<endl;
  if (frac>.5) frac = 1.-frac;
  int percent = (int)(frac * 100. + .5);
  float sdnorm = sd_baf_percent[50];
  float sdtumor = sd_baf_percent[percent];
  float var = pow((total_alpha)*sdnorm,2)
   + pow((1-total_alpha)*sdtumor,2);
  //cerr<<"debug: "<<pow((1-global_alpha-alpha)*sdtumor,2)<<endl;
  //cerr<<"alpha "<<global_alpha<<","<<alpha<<" SDs "<<sdnorm<<","<<sdtumor<<"; variance: "<<var<<endl;
  return sqrt(var);
}

float HMM_CNV::emission_prob(int state,float lrr, float baf,float pfb,ostream & os){
  float prob = 1; 
  float baf_noise = .01;
  float lrr_noise= .01;
  float alpha = state_total_alpha[state];
  float alpha_prob = state_alpha_prior[state];
  int tumorcn = state_tumorcn_vec[state];
  float lrr_prob = math->normal_pdf(lrr,state_mu_lrr[state]+sample_mu_lrr_adjusted[current_sample],state_sd_lrr[state]*sd_lrr_inflation);
  //cerr<<"Final SD and prob for state "<<state<<": "<<state_sd_lrr[state]*sd_lrr_inflation<<","<<lrr_prob<<endl;
  lrr_prob = lrr_noise+(1.-lrr_noise)*lrr_prob;

  float baf_prob = 0;
  float pointmass=.5;
  if (tumorcn==0){
    if (baf==0||baf==1){
      baf_prob+=normal_cdf;
    }else{
      baf_prob+=math->normal_pdf(baf,.5,sd_baf_homo_del*sd_baf_inflation);
    }
  }else{  // these other cases take into account stromal contamination
    bool status_loh = (tumorcn==1 || tumorcn==2 && state_het_ratio[state]==0);
    int effective_tumorcn = status_loh?1:tumorcn;
    // the first two cases only happen when tumor and normal are homozygote
    if (baf==0){
      baf_prob+=math->binom_prob(effective_tumorcn,0,pfb)*pointmass;
    }else if (baf==1){
      baf_prob+=math->binom_prob(effective_tumorcn,effective_tumorcn,pfb)
      * pointmass;
    }else{
      // this loop assumes the normal genotype is a heterozygote.
      for (int bac=0;bac<=effective_tumorcn;++bac){
        baf_prob+=math->binom_prob(effective_tumorcn,bac,pfb)*math->normal_pdf(baf,get_baf_mean(bac,effective_tumorcn,alpha),get_baf_sd(state,bac)*sd_baf_inflation);
      }
      if (status_loh){
        // assume homozygote normal too
        int bac = (baf<.5)?0:1;
        baf_prob+=math->binom_prob(effective_tumorcn,bac,pfb)*math->normal_pdf(baf,bac,sd_baf_percent[0]*sd_baf_inflation);
      }
    }
  }
  baf_prob = baf_noise+(1.-baf_noise)*baf_prob;
  //prob = lrr_prob ;
  prob = lrr_prob * baf_prob * alpha_prob;
  //prob = prob<HMM_CNV::small ?HMM_CNV::small:prob;
  bool debug = true;
  if (debug){
    os<<"ID:"<<state<<",";
    os<<"BAF:"<<baf<<",RR:"<<lrr<<",";
    os<<"P_RR:"<<lrr_prob<<",";
    os<<"P_BAF:"<<baf_prob<<",";
    os<<"P:"<<prob<<endl;
  }
  return prob;
}

state_t::state_t(int id,int tumor_cn, float global_alpha,float alpha, float het_ratio, bool alpha_identifiable,bool het_identifiable, float * sd_lrr_het, float * sd_lrr_loh, float  sd_lrr_normal){
  this->id = id;
  this->tumor_cn = tumor_cn;
  this->het_ident = het_identifiable;
  this->alpha_ident = alpha_identifiable;
  this->het_ratio = het_ratio;
  this->total_alpha = alpha;
  //this->total_alpha = (tumor_cn==0||tumor_cn==2) ? .01 : global_alpha+alpha*(1-global_alpha);
  float intercept = .32;
  float slope = (1.-intercept)/2.;
  this->mu_lrr = log2((1-total_alpha)*(intercept+tumor_cn*slope)+ (total_alpha));
  if (tumor_cn==2){
    if (het_ratio>0){
      // the variance of the normal tissue
      sd_lrr = (sd_lrr_normal);
    }else{
      sd_lrr = (sd_lrr_loh[tumor_cn]);
    }
  }else{
    sd_lrr = sd_lrr_het[tumor_cn];
  }
  //float var = pow(total_alpha*sd_lrr_normal,2);
  //cerr<<"id "<<id<<" first var "<<var<<" sdlrrnormal "<<sd_lrr_normal<<endl;
  //if (het_ratio>0){
    //var += pow((1-total_alpha)*sd_lrr_het[tumor_cn],2);
    //cerr<<"id "<<id<<" second  var "<<pow((1-total_alpha)*sd_lrr_het[tumor_cn],2)<<" sdlrrhet "<<sd_lrr_het[tumor_cn]<<endl;
  //}else{
    //var += pow((1-total_alpha)*sd_lrr_loh[tumor_cn],2);
  //}
  //cerr<<"id "<<id<<" final  var "<<var<<endl;
  //this->sd_lrr = sqrt(var);
}


void state_t::debug(){
  cerr<<id<<"\t"<<tumor_cn<<"\t"<<het_ratio<<"\t"<<total_alpha<<"\t"<<alpha_ident<<"\t"<<het_ident<<"\t"<<mu_lrr<<"\t"<<sd_lrr<<"\t"<<alpha_prior;
  cerr<<endl;
}

void HMM_CNV::init_state_objects(int current_sample){
  global_alpha = global_alphas[current_sample];
  int states = 0 ;
  float fraction = global_alpha;
  state_object_vector.clear();
  state_object_vector.push_back(state_t(states,0,global_alpha,fraction,HET_RATIO,false,false, sd_lrr_het,sd_lrr_loh,sd_lrr_normal));
  ++states;
  // normal state
  state_object_vector.push_back(state_t(states,2,global_alpha,fraction,HET_RATIO,false,true, sd_lrr_het,sd_lrr_loh,sd_lrr_normal));
  normal_state = states;
  ++states;
  // copy neutral LOH
  state_object_vector.push_back(state_t(states,2,global_alpha,fraction,0,false,true, sd_lrr_het,sd_lrr_loh,sd_lrr_normal));
  ++states;
  for (fraction=.01;fraction<1;fraction+=.25){
    // hemizygous deletion
    state_object_vector.push_back(state_t(states,1,global_alpha,fraction,HET_RATIO,true,false, sd_lrr_het,sd_lrr_loh,sd_lrr_normal));
    ++states;
    // amplifications
    for(int tumor_cn=3;tumor_cn<=4;++tumor_cn){
      state_object_vector.push_back(state_t(states,tumor_cn,global_alpha,fraction,HET_RATIO,true,false, sd_lrr_het,sd_lrr_loh,sd_lrr_normal));
      ++states;
    }
  }
  // States created with appropriate means and SDs. 
  cerr<<"state"<<"\t"<<"cn"<<"\t"<<"htzyg"<<"\t"<<"alpha"<<"\t"<<"alpha?\thet?\tmu_lrr\tsd_lrr\talpha_prior\n";
  for(int state=0;state<states;++state){
    state_mu_lrr[state] = state_object_vector[state].mu_lrr;
    state_het_ratio[state] = state_object_vector[state].het_ratio;
    state_total_alpha[state] = state_object_vector[state].total_alpha;
    state_tumorcn_vec[state] = state_object_vector[state].tumor_cn;
    state_het_ident[state] = state_object_vector[state].het_ident;
    state_alpha_ident[state] = state_object_vector[state].alpha_ident;
    state_sd_lrr[state] = state_object_vector[state].sd_lrr;
    state_alpha_prior[state] = state_object_vector[state].alpha_prior = 
    math->normal_pdf(state_total_alpha[state],global_alpha,sd_alpha); 
    state_object_vector[state].debug();
  }
}

void HMM_CNV::compute_emission_cache_mat(bool initgpu){
  float pfb = .5; // temporarily set this to .5
  //ofstream ofs_emission("debugging/emission_cache");
  ostringstream oss;
  oss<<"debugging/emission.debug."<<current_sample;
  ofstream ofs_debug(oss.str().data());
  for (int i=0;i<total_markers;++i){
    //if (i%1000==0) cerr<<".";
    bool debug = (iteration==-1)?true:false;
    float & lrr = lrr_vec[current_sample*total_markers+i];
    float & baf = baf_vec[current_sample*total_markers+i];
    if (debug) cerr<<"Observation: "<<i<<": "<<lrr<<","<<baf<<endl;
    float normalizer = 0;
    float raw_emission_prob[states];
    float max = 0;
    int maxindex;
    //float rr = pow(2,lrr);
    for (int j=0;j<states;++j){
      if (debug) cerr<<"State: "<<j<<endl;
      ofs_debug<<"Marker:"<<i<<",";
      raw_emission_prob[j] = emission_prob(j,lrr,baf,pfb,ofs_debug);
      if (debug) cerr<<"Raw prob: "<<raw_emission_prob[j]<<endl<<endl;
      if (initgpu){
        //raw_emission_prob[j] = emission_rr[i*states+j];
      }
      //raw_emission_prob[j] = state_object_vector[j].emission_prob(lrr_vec[i],baf_vec[i],sd_rr,sd_baf,alpha_weight_map,alpha_block_vec[i],debug);
      //if(j) ofs_emission<<"\t";
      //ofs_emission<<raw_emission_prob[j];
      if (raw_emission_prob[j]>max) {
        max = raw_emission_prob[j];
        maxindex = j;
      }
      normalizer+=raw_emission_prob[j];
    }
    //state_t s = state_object_vector[maxindex];
    //cerr<<"Max at state "<<maxindex<<" fractional: "<<s.mixture_cn<<" normal BAF "<<1.*s.normal_bac/s.normal_cn<<" tumor BAF "<<1.*s.tumor_bac/s.tumor_cn<<endl;
    for (int j=0;j<states;++j){
      float e = raw_emission_prob[j]/normalizer;
      if (e<small) e = small;
      emission_cache_mat[current_sample*hmm_matrix_size+i*states+j] = e;
    }
    //ofs_emission<<endl;
  }
  //ofs_emission.close();
  ofs_debug.close();
  for(int j=0;j<states;++j){
    forward_mat_log[current_sample*hmm_matrix_size+j] = log(emission_cache_mat[current_sample*hmm_matrix_size+j]);
    backward_mat_log[current_sample*hmm_matrix_size+(total_markers-1)*states+j] = 0;
  }
  cerr<<"Done initializing emission_cache_mat\n";
  if (use_gpu && initgpu){
     #ifdef USE_GPU
     cerr<<"Transferring matrices to GPU\n";
     err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), forward_mat_log , NULL, NULL );
    clSafe(err, "Transfer forward mat to GPU");
     err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), backward_mat_log , NULL, NULL );
    clSafe(err, "Transfer backward mat to GPU");
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.emission_cache_mat_buffer, CL_TRUE, 0,  sizeof(float)*(total_markers*states), emission_cache_mat , NULL, NULL );
    clSafe(err, "Transfer emission cache mat to GPU");
    //err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.greedy_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), greedy_mat_log , NULL, NULL );
    //clSafe(err, "Transferring greedy  mat to GPU");
     cerr<<"Transferred matrices to GPU\n";
    #endif
  }
}


double HMM_CNV::addMultipleLogs(double * logs, int len){
  //cerr<<"Called addMultipleLogs\n";
  double avg=0;
  for(int i=0;i<len;++i){
    avg-=logs[i];
  }
  avg/=len;
  //avg-=10;
  //cerr<<"Average scaling: "<<avg<<endl;
  double val = 0;
  for(int i=0;i<len;++i){
    val+=exp(logs[i]+avg);
  }
  double a = log(val)-avg;
  //cerr<<"Results "<<a<<endl;
  return a;
}

//float HMM_CNV::addTwoLogs(float log1,float log2){
//  float logarr[2];
//  logarr[0] = log1;
//  logarr[1] = log2;
//  return addMultipleLogs(logarr,2);
//}

void HMM_CNV::forwardbackward_gpu(){
  #ifdef USE_GPU
  int counter[states];
  for(int i=0;i<states;++i) counter[i] = 1;
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.counter_index_buffer, CL_TRUE, 0,  states*sizeof(int), &counter , NULL, NULL );
  clSafe(err, "writing buffer for counter index");
  cerr<<"ForwardBackward\n";
  for(int obs=1;obs<total_markers;++obs){
    float forwardgap = gap_vec[obs];
    float backwardgap = gap_vec[total_markers-obs];
    //cerr<<"1 Writing for obs "<<obs<<endl;
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_gap_buffer, CL_TRUE, 0,  sizeof(float), &forwardgap , NULL, NULL );
    clSafe(err,"write forward gap");
    //cerr<<"2 Writing for obs "<<obs<<endl;
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_gap_buffer, CL_TRUE, 0,  sizeof(float), &backwardgap , NULL, NULL );
    clSafe(err,"write backward gap");
    cl::Event event_rescale;
    //cerr<<"1 Launching for obs "<<obs<<endl;
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_rescale_transition,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,2),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&event_rescale);
    event_rescale.wait();
    clSafe(err,"CommandQueue::enqueueNDRangeRescale Transition()");
    //obs=-1;
    bool debug_scaled = false;
    if (debug_scaled){
      float test[states*states];
      err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.scaled_forward_transition_mat_buffer, CL_TRUE, 0, (states*states)*sizeof(float),test,NULL,NULL);
      clSafe(err,"CommandQueue::Read scaled forward()");
      ostringstream oss;
      oss<<"debugging/scaled.gpu."<<obs;
      ofstream ofs_scaled(oss.str().data());
      for(int i=0;i<states;++i){
        for(int j=0;j<states;++j){
          if (j) ofs_scaled<<" ";
          ofs_scaled<<j<<":"<<test[i*states+j];
        }
        ofs_scaled<<endl;
      }
      ofs_scaled.close();
      exit(0);
    }

    //cerr<<"1 Writing for obs "<<obs<<endl;
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_gap_buffer, CL_TRUE, 0,  sizeof(float), &forwardgap , NULL, NULL );
    clSafe(err,"write forward gap");
    cl::Event event_forward_backward;
  
    //cerr<<"2 Launching for obs "<<obs<<endl;
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_forward_backward,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,2),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&event_forward_backward);
    clSafe(err,"CommandQueue::enqueueNDRangeKernelForwardBackward()");
    event_forward_backward.wait();
    //cerr<<"Done\n";
    if (obs % 1000 == 0) cerr<<" "<<obs;
  }
  cerr<<"Launched forward backward\n";
  cerr<<endl;
  err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_likelihood,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,NULL);
  clSafe(err,"CommandQueue::enqueueNDRangeKernellikelihood()");
  float likelihood_log;
  err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.likelihood_log_buffer, CL_TRUE, 0, (1)*sizeof(float),&likelihood_log,NULL,NULL);
  clSafe(err,"CommandQueue::Readlikelihood()");
  cerr<<"Iteration: "<<iteration<<": GPU likelihood_log: "<<likelihood_log<<endl<<endl;
  bool debug = false;
  if (debug ){
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.forward_mat_buffer, CL_TRUE, 0, (total_markers*states)*sizeof(double),forward_mat_log);
    clSafe(err,"CommandQueue::readforward()");
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backward_mat_buffer, CL_TRUE, 0, (total_markers*states)*sizeof(double),backward_mat_log);
    clSafe(err,"CommandQueue::read prob Backward()");
    cerr<<"GPU forward/backward answer:\n";
    for(int i=0;i<total_markers;++i){
      cerr<<"obs: "<<i;
      for(int j=0;j<states;++j){
        cerr<<" state:"<<j<<": "<<forward_mat_log[i*states+j]<<"/"<<backward_mat_log[i*states+j]<<" "<<forward_mat_log[i*states+j]+backward_mat_log[i*states+j]-likelihood_log<<endl;
      }
      cerr<<endl;
    }
    exit(0);
  }
  // COMPUTE PROBABILITY MATRIX
  err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_posterior,cl::NullRange,cl::NDRange(total_markers*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,NULL);
  //eventlist.push_back(e);
  clSafe(err,"CommandQueue::enqueueNDRangeKernelPosterior()");
  bool debug2 = false;
  if (debug2){
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.prob_mat_buffer, CL_TRUE, 0, (total_markers*states)*sizeof(double),prob_mat_log,NULL,NULL);
    clSafe(err,"CommandQueue::readbuffer prob mat");
    cerr<<"\nGPU: prob mat\n";
    for(int i=0;i<total_markers;++i){
       cerr<<"obs "<<i<<":";
       for(int j=0;j<states;++j){
        //cerr<<" "<<i<<":"<<prob_mat_log[j*total_markers+i]<<endl;
        cerr<<" "<<prob_mat_log[i*states+j];
      }
      cerr<<endl<<endl;
    }
    exit(0);
  }
  #endif
}

void HMM_CNV::marginalize_posteriors(){
  ostringstream oss1;
  ostringstream oss2;
  ostringstream oss3;
  ostringstream oss4;
  ostringstream oss5;
  oss1<<"results/sample."<<sample_annotations[current_sample];
  oss2<<"debugging/posterior.alpha."<<current_sample;
  oss3<<"debugging/posterior.het."<<current_sample;
  oss4<<"debugging/posterior.tumor_cn."<<current_sample;
  oss5<<"debugging/transition."<<current_sample;
  cerr<<"Marginalizing for sample "<<current_sample<<" with "<<total_markers<<" markers"<<endl;
  ofstream ofs_expected(oss1.str().data());
  ofs_expected<<"chr\tpos\tBAF\tRR\ttumor_cn\talpha\theterozyg\n";
  ofstream ofs_alpha(oss2.str().data());
  ofstream ofs_het(oss3.str().data());
  ofstream ofs_tumor_cn(oss4.str().data());
  ofstream ofs_trans(oss5.str().data());
  for(int i=0;i<states;++i){
    for(int j=0;j<states;++j){
      if (j) ofs_trans<<"\t";
      ofs_trans<<sample_trans_mat[current_sample*trans_matrix_size+i*states+j];
    }
    ofs_trans<<endl;
  }
  ofs_trans.close();
  float alpha_probs[100];
  float het_probs[100];
  int tumor_cn_elements = 5;
  float tumor_cn_probs[tumor_cn_elements];
  for(int i = 0;i<total_markers;++i){
    memset(alpha_probs,0,sizeof(alpha_probs));
    memset(het_probs,0,sizeof(het_probs));
    memset(tumor_cn_probs,0,sizeof(tumor_cn_probs));
    for(int state=0;state<states;++state){
      state_t s = state_object_vector[state];
      float prob = exp(prob_mat_log[current_sample*hmm_matrix_size+i*states+state]);
      tumor_cn_probs[s.tumor_cn]+=prob;
      if(s.alpha_ident){
        int alpha_int = (int)(s.total_alpha*100);
        alpha_probs[alpha_int]+=prob;
      }
      if(s.het_ident){
        int het_int = (int)(s.het_ratio*100);
        het_probs[het_int]+=prob;
      }
    }
    float alpha_normalizer = 0;
    float het_normalizer = 0;
    for(int j=0;j<100;++j){
      if (alpha_probs[j]>0) alpha_normalizer += alpha_probs[j];
      if (het_probs[j]>0) het_normalizer += het_probs[j];
    }
    ofs_alpha<<i;
    ofs_het<<i;
    float exp_alpha = 0;
    float exp_het = 0;
    for(int j=0;j<100;++j){
      if (alpha_probs[j]>0){
         alpha_probs[j]/=alpha_normalizer;
         ofs_alpha<<"\t"<<(j/100.)<<":"<<alpha_probs[j];
         exp_alpha+=alpha_probs[j] * (j/100.);
      }
      if (het_probs[j]>0){
         het_probs[j]/=het_normalizer;
         ofs_het<<"\t"<<(j/100.)<<":"<<het_probs[j];
         exp_het+=het_probs[j] * (j/100.);
      }
    }
    ofs_alpha<<endl;
    ofs_het<<endl;
    float normalizer = 0;
    for(int j=0;j<tumor_cn_elements;++j){
      if (tumor_cn_probs[j]>0) normalizer += tumor_cn_probs[j];
    }
    ofs_tumor_cn<<i;
    float exp_tumor_cn=0;
    for(int j=0;j<tumor_cn_elements;++j){
      if (tumor_cn_probs[j]>0) {
        tumor_cn_probs[j]/=normalizer;
        ofs_tumor_cn<<"\t"<<j<<":"<<tumor_cn_probs[j];
        exp_tumor_cn+=tumor_cn_probs[j] * j;
      }
    }
    ofs_tumor_cn<<endl;
    ofs_expected<<chr<<"\t"<<pos_vec[i]<<"\t"<<baf_vec[current_sample*total_markers+i]<<"\t"<<lrr_vec[current_sample*total_markers+i]<<"\t"<<exp_tumor_cn<<"\t"<<exp_alpha<<"\t"<<exp_het<<endl;
  }
  ofs_tumor_cn.close();
  ofs_alpha.close();
  ofs_het.close();
  ofs_expected.close();
}

//void HMM_CNV::viterbi_output(){
//  
//  bool output_gpu = true;
//  bool output_penncnv = false;
//  bool output_verbose = false;
//  ofstream ofs1;
//  if (output_gpu){
//    int lastid = -1;
//    int startindex = 0;
//    int endindex = -1;
//    ofs1.open("results.gpu");
//    ofs1<<"chr\tstart\tend\tlength\tstate\tcn\tbac\tsnp1\tsnp2\talpha\n";
//    for(int i = 0;i<total_markers;++i){
//      state_t s = state_object_vector[backtrace_vec[i]];
//      if (i && s.id!=lastid){
//        if (endindex>startindex){
//          state_t s2 = state_object_vector[backtrace_vec[endindex]];
//          float a = s2.tumor_cn==-2?0.5:s2.alpha;
//          ofs1<<chr<<"\t"<<pos_vec[startindex]<<"\t"<<pos_vec[endindex]<<"\t"<<1.*(pos_vec[endindex]-pos_vec[startindex])/1e6<<"\t"<<s2.id<<"\t"<<s2.tumor_cn<<"\t"<<s2.tumor_bac<<"\t"<<rslist[startindex]<<"\t"<<rslist[endindex]<<"\t"<<a<<endl;
//        }
//        startindex = i;
//      }else{
//        endindex = i;
//      }
//      lastid = s.id;
//    }
//    state_t s2 = state_object_vector[backtrace_vec[endindex]];
//    float a = s2.tumor_cn==-2?0.5:s2.alpha;
//    ofs1<<chr<<"\t"<<pos_vec[startindex]<<"\t"<<pos_vec[endindex]<<"\t"<<1.*(pos_vec[endindex]-pos_vec[startindex])/1e6<<"\t"<<s2.id<<"\t"<<s2.tumor_cn<<"\t"<<s2.tumor_bac<<"\t"<<rslist[startindex]<<"\t"<<rslist[endindex]<<"\t"<<a<<endl;
//    ofs1.close();
//  }
//  if (output_penncnv){
//    ofs1.open("results.penncnv");
//    int lastid = -1;
//    int startindex = 0;
//    int endindex = -1;
//    for(int i = 0;i<total_markers;++i){
//      state_t s = state_object_vector[backtrace_vec[i]];
//      if (i && s.id!=lastid){
//        if (endindex>startindex){
//          state_t s2 = state_object_vector[backtrace_vec[endindex]];
//          int numsnp = endindex-startindex+1;
//          int length = pos_vec[endindex]-pos_vec[startindex]+1;
//          ofs1<<"chr"<<chr<<":"<<pos_vec[startindex]<<"-"<<pos_vec[endindex]<<"\tnumsnp="<<numsnp<<"\tlength="<<length<<"\tstate="<<s2.id<<",cn="<<s2.tumor_cn<<",bac="<<s2.tumor_bac<<"\talpha="<<s2.alpha<<"\tstartsnp="<<rslist[startindex]<<"\tendsnp="<<rslist[endindex]<<"\t0\t0\n";
//        }
//        startindex = i;
//      }else{
//        endindex = i;
//      }
//      lastid = s.id;
//    }
//    state_t s2 = state_object_vector[backtrace_vec[endindex]];
//    int numsnp = endindex-startindex+1;
//    int length = pos_vec[endindex]-pos_vec[startindex]+1;
//    ofs1<<"chr"<<chr<<":"<<pos_vec[startindex]<<"-"<<pos_vec[endindex]<<"\tnumsnp="<<numsnp<<"\tlength="<<length<<"\tstate="<<s2.id<<",cn="<<s2.tumor_cn<<",bac="<<s2.tumor_bac<<"\talpha="<<s2.alpha<<"\tstartsnp="<<rslist[startindex]<<"\tendsnp="<<rslist[endindex]<<"\t0\t0\n";
//    ofs1.close();
//  }
//  if (output_verbose){
//    ostringstream oss_1;
//    oss_1<<"results.verbose.iter."<<iteration;
//    ofs1.open(oss_1.str().data());
//    ofs1<<"index\tsnp\tchr\tposition\tstate_id\talpha\tbac_norm\tcn_tumor\tbac_tumor\n";
//    //ofs1<<"index\tsnp\tchr\tposition\tstate_id\talpha\tbac_norm\tcn_tumor\tbac_tumor\temission_baf\temission_rr\tnormal_rr\t0\t.5\t1\n";
//    for(int i = 0;i<total_markers;++i){
//      state_t s = state_object_vector[backtrace_vec[i]];
//      float alpha = s.alpha;
//      ofs1<<i<<"\t"<<rslist[i]<<"\t"<<chr<<"\t"<<pos_vec[i]<<"\t"<<s.id<<"\t"<<alpha<<"\t"<<s.normal_bac<<"\t"<<s.tumor_cn<<"\t"<<s.tumor_bac<<endl;
//      //ofs1<<i<<"\t"<<rslist[i]<<"\t"<<chr<<"\t"<<pos_vec[i]<<"\t"<<s.id<<"\t"<<alpha<<"\t"<<s.normal_bac<<"\t"<<s.tumor_cn<<"\t"<<s.tumor_bac<<"\t"<<emission_baf[i*states+backtrace_vec[i]] <<"\t"<<emission_rr[i*states+backtrace_vec[i]]<<"\t"<<emission_rr[i*states]<<"\t"<<emission_baf[i*states]<<"\t"<<emission_baf[i*states+9]<<"\t"<<emission_baf[i*states+18]<<endl;
//    }
//    ofs1.close();
//  }
//}


//void HMM_CNV::viterbi_gpu(){
//  #ifdef USE_GPU
//  int counter[states];
//  for(int i=0;i<states;++i) counter[i] = 1;
//  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.counter_index_buffer, CL_TRUE, 0,  states*sizeof(int), &counter , NULL, NULL );
//  clSafe(err, "writing buffer for counter index");
//
//  double LOG_MIN = -1e11;
//  double max = LOG_MIN;
//  int obs = 0;
//  for(int j=0;j<states;++j){
//     //cerr<<" "<<greedy_mat_log[j];
//     if (greedy_mat_log[j] > max){
//       backtrace_vec[obs] = j;
//       max = greedy_mat_log[j];
//     }
//  }
//  //exit(0);
//  //err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backtrace_vec_buffer, CL_TRUE, 0,  total_markers*sizeof(int), &backtrace_vec , NULL, NULL );
//  //clSafe(err, "writing buffer for backtrace index");
//
//  for(int obs=1;obs<total_markers;++obs){
//    float forwardgap = gap_vec[obs];
//    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_gap_buffer, CL_TRUE, 0,  sizeof(float), &forwardgap , NULL, NULL );
//    clSafe(err,"CommandQueue::update forward gap()");
//    cl::Event rescaleEvent,viterbiEvent;
//    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_rescale_transition,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&rescaleEvent);
//    clSafe(err,"CommandQueue::enqueueNDRangeRescale Transition()");
//    rescaleEvent.wait();
//    //cerr<<"Rescaled for obs: "<<obs<<endl;
//    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_viterbi,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&viterbiEvent);
//    clSafe(err,"CommandQueue::enqueueNDRangeViterbi()");
//    viterbiEvent.wait();
//    //cerr<<"Viterbi for obs: "<<obs<<endl;
//  }
//  err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.greedy_mat_buffer, CL_TRUE, 0, (total_markers*states)*sizeof(double),greedy_mat_log,NULL,NULL);
//  clSafe(err,"CommandQueue::readbuffer greedy mat");
//  err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.bestpath_mat_buffer, CL_TRUE, 0, (total_markers*states)*sizeof(int),bestpath_mat,NULL,NULL);
//  clSafe(err,"CommandQueue::readbuffer bestpath mat");
//  err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.backtrace_vec_buffer, CL_TRUE, 0, (total_markers)*sizeof(int),backtrace_vec,NULL,NULL);
//  clSafe(err,"CommandQueue::readbuffer backtrace vec");
//  backtrace_vec[total_markers-1] = 0;
//  double currval = LOG_MIN;
//  for(int state=0;state<states;++state){
//    if (greedy_mat_log[(total_markers-1)*states+state] > currval){
//      currval = greedy_mat_log[(total_markers-1)*states+state];
//      backtrace_vec[total_markers-1] = state;
//    }
//  }
//  for(int obs = total_markers-2;obs>=0;--obs){
//    backtrace_vec[obs] = bestpath_mat[(obs+1)*states+backtrace_vec[obs+1]];
//  }
//  viterbi_output();
//  #endif
//}

//void HMM_CNV::viterbi_cpu(){
//  if (use_gpu){
//    #ifdef USE_GPU
//    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.transition_mat_buffer, CL_TRUE, 0, (states*states)*sizeof(float),trans_mat);
//    clSafe(err,"CommandQueue::readforward()");
//    #endif
//  }
//  //float sum,forwardsum,backwardsum;
//  // transpose the transition matrix to avoid strided memory access
//  //float trans_mat_t[states*states];
//  for(int i=0;i<states;++i){
//    for(int j=0;j<states;++j){
//      //trans_mat_t[j*states+i] = trans_mat[i*states+j];
//    }
//  }
//  //int back_trace[total_markers][states];
//  //int trace[total_markers];
//  double LOG_MIN = -1e11;
//  double max = LOG_MIN;
//  int obs = 0;
//  for(int j=0;j<states;++j){
//     //cerr<<" "<<greedy_mat_log[j];
//     if (greedy_mat_log[j] > max){
//       backtrace_vec[obs] = j;
//       max = greedy_mat_log[j];
//     }
//  }
//  //cerr<<"Max is "<<max<<endl;
//  for(int obs=1;obs<total_markers;++obs){
//    float new_trans_mat[states*states];
//    scaleTransitionMat(new_trans_mat,gap_vec[obs]);
//    for(int currentstate=0;currentstate<states;++currentstate){
//      double max = LOG_MIN;
//      int maxind = 0;
//      for(int prevstate=0;prevstate<states;++prevstate){
//    //cerr<<"prevstate "<<prevstate<<" currstate "<<currentstate<<" Got here\n";
//        //double curr_val = greedy_mat_log[(obs-1)*states+prevstate] ;
//        double curr_val = greedy_mat_log[(obs-1)*states+prevstate] + log(new_trans_mat[prevstate*states+currentstate]);
//    //cerr<<"prevstate "<<prevstate<<" currstate "<<currentstate<<" Got here\n";
//        if (curr_val>=max){
//          maxind = prevstate;
//          max = curr_val;
//        }
//      }
//      greedy_mat_log[obs*states+currentstate] = max + log(emission_cache_mat[obs*states+currentstate]);
//      //cerr<<obs<<" "<<currentstate<<" Max state "<<maxind<<" max prob: "<<max<<" greedy "<<greedy_mat_log[obs*states+currentstate]  <<" emission: "<<log(emission_cache_mat[obs*states+currentstate]) <<endl;
//      bestpath_mat[obs*states+currentstate] = maxind;
//    }
//  }
//  backtrace_vec[total_markers-1] = 0;
//  double currval = LOG_MIN;
//  for(int state=0;state<states;++state){
//    if (greedy_mat_log[(total_markers-1)*states+state] > currval){
//      backtrace_vec[total_markers-1] = state;
//      currval = greedy_mat_log[(total_markers-1)*states+state];
//    }
//  }
//  for(int obs = total_markers-2;obs>=0;--obs){
//    backtrace_vec[obs] = bestpath_mat[(obs+1)*states+backtrace_vec[obs+1]];
//  }
//  viterbi_output();
//}

void HMM_CNV::forwardbackward_cpu(){
  //ofstream ofs_forward("debugging/forward");
  // RUN FORWARD ALGORITHM
  for(int obs=1;obs<total_markers;++obs){
  //cerr<<"Launching obs"<<total_markers<<"\n";
    if (obs % 1000 == 0) cerr<<" "<<obs;
    float new_trans_mat[states*states];
    scaleTransitionMat(new_trans_mat,gap_vec[obs]);
    for(int currentstate=0;currentstate<states;++currentstate){
      double sum_log[states];
      for(int prevstate=0;prevstate<states;++prevstate){
        // OPTIMIZED
        //sum_log[prevstate]=  log(emission_cache_mat[obs*states+currentstate]);
        sum_log[prevstate]=forward_mat_log[current_sample*hmm_matrix_size+(obs-1)*states+prevstate] + log(new_trans_mat[prevstate*states+currentstate]) + log(emission_cache_mat[current_sample*hmm_matrix_size+obs*states+currentstate]);
      }
      forward_mat_log[current_sample*hmm_matrix_size+obs*states+currentstate] = addMultipleLogs(sum_log,states);
      //if (currentstate) ofs_forward<<"\t";
      //ofs_forward<<forward_mat_log[current_sample*hmm_matrix_size+obs*states+currentstate];
      //if (obs>-1 && obs<102) cerr<<"obs: "<<obs<<" cur: "<<currentstate<<" new prob: "<<forward_mat_log[(obs)*states+currentstate]<<endl;
    }
    //ofs_forward<<endl;
  }
  //ofs_forward.close();
  //exit(0);
  // initialize the probabilities at the last observation as 1 since no data beyond
  // RUN BACKWARD ALGORITHM
  for(int obs=total_markers-2;obs>=0;--obs){
    if (obs % 1000 == 0) cerr<<" "<<obs;
    float new_trans_mat[states*states];
    scaleTransitionMat(new_trans_mat,gap_vec[obs+1]);
    for(int currentstate=0;currentstate<states;++currentstate){
      double sum_log[states];
      for(int futurestate=0;futurestate<states;++futurestate){
        // OPTIMIZED
        //sum_log[futurestate] = log(emission_cache_mat[(obs+1)*states+futurestate]); 
        sum_log[futurestate] = log(new_trans_mat[currentstate*states+futurestate]) +log(emission_cache_mat[current_sample*hmm_matrix_size+(obs+1)*states+futurestate]) + backward_mat_log[current_sample*hmm_matrix_size+(obs+1)*states+futurestate];
      }
      backward_mat_log[current_sample*hmm_matrix_size+obs*states+currentstate] = addMultipleLogs(sum_log,states);
    }
  }
  double sum_log[states];
  for(int state=0;state<states;++state){
    sum_log[state]=forward_mat_log[current_sample*hmm_matrix_size+state]+backward_mat_log[current_sample*hmm_matrix_size+state];
  }
  float likelihood_log = addMultipleLogs(sum_log,states);
  cerr<<"Iteration: "<<iteration<<": CPU likelihood_log: "<<likelihood_log<<endl<<endl;
  sample_likelihood_log[current_sample] = likelihood_log;
  //exit(0);
  bool debug = false;
  if (debug){
    for(int i=0;i<total_markers;++i){
      cerr<<"obs: "<<i<<endl;
      for(int j=0;j<states;++j){
        cerr<<"state: "<<j<<": "<<forward_mat_log[current_sample*hmm_matrix_size+i*states+j]<<"/"<<backward_mat_log[current_sample*hmm_matrix_size+i*states+j]<<" "<<forward_mat_log[current_sample*hmm_matrix_size+i*states+j]+backward_mat_log[current_sample*hmm_matrix_size+i*states+j]-likelihood_log<<endl;
      }
      cerr<<endl;
    }
    exit(0);
  }
  for(int obs=0;obs<total_markers;++obs){
    for(int state=0;state<states;++state){
      //prob_mat_log_t[state*total_markers+obs] = forward_mat_log_t[state*total_markers+obs] +  backward_mat_log_t[state*total_markers+obs]-likelihood_log; 
      prob_mat_log[current_sample*hmm_matrix_size+obs*states+state] = forward_mat_log[current_sample*hmm_matrix_size+obs*states+state] +  backward_mat_log[current_sample*hmm_matrix_size+obs*states+state]-likelihood_log; 
    }
  }
  bool debug2 = false;
  if (debug2){
    cerr<<"CPU Probability matrix: \n";
    for(int obs=0;obs<total_markers;++obs){
      //float rowsum=0;
      cerr<<"obs:"<<obs;
      for(int state=0;state<states;++state){
        cerr<<" "<<prob_mat_log[current_sample*hmm_matrix_size+obs*states+state];
        //rowsum+=exp(prob_mat_log[obs*states+state]);
      }
    cerr<<endl;
    //cerr<<" "<<rowsum<<endl;
    }
    exit(0);
  }
}

void HMM_CNV::baumwelch_gpu(){
  #ifdef USE_GPU
  cl_int err;
  // COMPUTE THE NORMALIZING FACTOR
  err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_baumwelch_normalize,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_OBS_BLOCK_WIDTH,1),cl::NDRange(HMM_OBS_BLOCK_WIDTH,1),NULL,NULL);
  bool debug1 = false;
  if (debug1){
    double denom_f[states];
    for (int i=0;i<states;++i) denom_f[i]=0;
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.baumwelch_norm_buffer, CL_TRUE, 0, (states)*sizeof(double),denom_f);
    clSafe(err, "read norm vec");
    cerr<<"GPU denom:\n";
    for (int i=0;i<states;++i) {
      cerr<<" "<<denom_f[i];
    }
    cerr<<endl;
    exit(0);
  }
  // EXECUTE THE TRANSITION MATRIX UPDATE
  //clSafe(err," forward transition");
  cl::Event zeroEvent;
  err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_zero_transition,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&zeroEvent);
  clSafe(err, "kernel zero transition");
  zeroEvent.wait();
  int counter[states];
  for(int i=0;i<states;++i) counter[i] = 0;
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.counter_index_buffer, CL_TRUE, 0,  states*sizeof(int), &counter , NULL, NULL );
  for(int obs=0;obs<total_markers-1;++obs){
    float forwardgap = gap_vec[obs+1];
    err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_gap_buffer, CL_TRUE, 0,  sizeof(float), &forwardgap , NULL, NULL );
    clSafe(err,"CommandQueue::update forward gap()");
    cl::Event rescaleEvent,baumwelchEvent;
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_rescale_transition,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&rescaleEvent);
    clSafe(err,"CommandQueue::enqueueNDRangeRescale Transition()");
    rescaleEvent.wait();
    err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_baumwelch_transition,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,&baumwelchEvent);
    clSafe(err, "kernel transition");
    baumwelchEvent.wait();
    if (obs % 1000 == 0) cerr<<" "<<obs;
  }
  cerr<<endl;
  err = opencl_info.command_queue.enqueueNDRangeKernel(opencl_info.kernel_normalize_transition,cl::NullRange,cl::NDRange(HMM_STATE_BLOCK_WIDTH*HMM_STATE_BLOCK_WIDTH,1),cl::NDRange(HMM_STATE_BLOCK_WIDTH,1),NULL,NULL);
  clSafe(err, "kernel norm transition");
  bool debug2 = false;
  if (debug2){
    err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.transition_mat_buffer, CL_TRUE, 0, (states*states)*sizeof(float),trans_mat);
    clSafe(err, "read new transition");
    cerr<<"GPU transition\n";
    for(int i=0;i<states;++i){
      float rowsum=0;
      for(int j=0;j<states;++j){
        rowsum+=trans_mat[i*states+j];
        cerr<<" "<<trans_mat[i*states+j];
      }
      cerr<<" rowsum:"<<rowsum<<endl;
    }
    exit(0);
  }
  #endif
  // get MLE of SD
}

void HMM_CNV::baumwelch_cpu(){
  cerr<<"Running Baum Welch\n";
  //cerr<<"Old transition\n";
  for(int i=0;i<states;++i){
    float rowsum = 0;
    for(int j=0;j<states;++j){
      rowsum+=trans_mat[i*states+j];
      //cerr<<" "<<trans_mat[i*states+j];
    }
    //cerr<<" rowsum: "<<rowsum<<endl;

  }

  double denom_f_arr[states];
  for(int i=0;i<states;++i){
    denom_f_arr[i] = 0;
    //float test2 = 0;
    for(int obs=0;obs<total_markers-1;++obs){
      // OPTIMIZED
      denom_f_arr[i]+=exp(prob_mat_log[current_sample*hmm_matrix_size+obs*states+i]) ;
      //denom_f_arr[i]+=exp(prob_mat_log_t[i*total_markers+obs]) ;
    }
  }
  bool debug1 = false;
  if (debug1){
    cerr<<"CPU denom:\n";
    for(int i=0;i<states;++i){
      cerr<<" "<<denom_f_arr[i];
    }
    cerr<<endl;
    exit(0);
  }
  // new transition kernel
  double numer_f[states*states];
  for(int i=0;i<states*states;++i) numer_f[i]=0.;
  float likelihood_log=sample_likelihood_log[current_sample];
  for(int obs=0;obs<total_markers-1;++obs){
    if (obs % 1000 == 0) cerr<<" "<<obs;
    float new_trans_mat[states*states];
    scaleTransitionMat(new_trans_mat,gap_vec[obs+1]);
    for(int i=0;i<states;++i){
      for(int j=0;j<states;++j){
        //cerr<< "debug: "<<forward_mat_log[obs*states+i]<<","<<log(new_trans_mat[i*states+j])<<","<<log(emission_cache_mat[(obs+1)*states+j]) <<","<< backward_mat_log[(obs+1) * states + j] <<","<< likelihood_log<<endl;
        numer_f[i*states+j]+=exp(forward_mat_log[current_sample*hmm_matrix_size+obs*states+i]+log(new_trans_mat[i*states+j])+log(emission_cache_mat[current_sample*hmm_matrix_size+(obs+1)*states+j]) + backward_mat_log[current_sample*hmm_matrix_size+(obs+1) * states + j] - likelihood_log );
        //numer_f[i*states+j]+=exp(forward_mat_log_t[i*total_markers+obs]+log(new_trans_mat[i*states+j])+log(emission_cache_mat[(obs+1)*states+j]) + backward_mat_log_t[j*total_markers+(obs+1)] - likelihood_log );
      }
    }
    if (obs==-1){
    for(int i=0;i<states;++i){
      float rowsum = 0;
      for(int j=0;j<states;++j){
        rowsum+=numer_f[i*states+j];
        cerr<<" "<<numer_f[i*states+j];
      }
      cerr<<" rowsum: "<<rowsum<<endl;
    }
    exit(0);
    }
  }
  /**
    for(int i=0;i<states;++i){
      float rowsum = 0;
      for(int j=0;j<states;++j){
        rowsum+=numer_f[i*states+j];
        cerr<<" "<<numer_f[i*states+j];
      }
      cerr<<" rowsum: "<<rowsum<<endl;
    }
    exit(0);
  **/
  for(int i=0;i<states;++i){
    for(int j=0;j<states;++j){
      //cerr<<"numer: "<<numer_f[i*states+j]<<endl;
      trans_mat[i*states+j]=(numer_f[i*states+j]/denom_f_arr[i]);
      //trans_mat[i*states+j]= HMM_SMALL + (1.-HMM_SMALL)*trans_mat[i*states+j];
    }
  }
  
  bool debug2 = false;
  if (debug2){
    cerr<<"CPU transition\n";
    for(int i=0;i<states;++i){
      float rowsum = 0;
      for(int j=0;j<states;++j){
        rowsum+=trans_mat[i*states+j];
        cerr<<" "<<trans_mat[i*states+j];
      }
      cerr<<" rowsum: "<<rowsum<<endl;
    }
    exit(0);
  }
  // FOR EMISSIONS COMPUTE MLE SIGMAS
  
  
  
}

//void HMM_CNV::get_sd_mle(bool weighted){
//  //sd_baf = .47;
//  //sd_rr = .42
//  //sd_rr = sd_baf = .3;
//  //return;
//  if (use_gpu){
//    //err = opencl_info.command_queue.enqueueReadBuffer(opencl_info.prob_mat_buffer, CL_TRUE, 0, (total_markers*states)*sizeof(double),prob_mat_log,NULL,NULL);
//    //clSafe(err,"CommandQueue::readbuffer prob mat");
//  }
//  
//  float rr_numer = 0;
//  float baf_numer = 0;
//  float rr_denom = 0;
//  float denom_vec[states];
//  for(int i=0;i<states;++i){
//    denom_vec[i] = 0;
//    sd_rr_vec[i] = 0;
//    sd_baf_vec[i] = 0;
//  }
//  float baf_denom = 0;
//  for(int i=0;i<total_markers;++i){
//    float r = pow(2,lrr_vec[i]);
//    //for(int j=1;j<=1;++j){
//    for(int j=0;j<states;++j){
//      state_t s = state_object_vector[j];
//      float p = weighted?emission_cache_mat[i*states+j]:1;
//      //float p = weighted?exp(prob_mat_log[i*states+j]):1;
//      bool use_global = true;
//      float r_mean=0,baf_mean=0;
//      if(use_global){
//        r_mean = rr_offset;
//        baf_mean = .5;
//      }else{
//        r_mean = s.mu_rr;
//        //baf_mean = s.mu_baf;
//      }
//      rr_numer+=p*pow(r-r_mean,2);
//      baf_numer+=p*pow(baf_vec[i]-baf_mean,2);
//      rr_denom+=p;
//      baf_denom+=p;
//
//      sd_rr_vec[j]+=p*pow(r-r_mean,2);
//      sd_baf_vec[j]+=p*pow(baf_vec[i]-baf_mean,2);
//      denom_vec[j]+=p;
//    }
//  }
//  sd_baf = sqrt(baf_numer/baf_denom);
//  sd_rr = sqrt(rr_numer/rr_denom);
//  cerr<<"MLE SDs for global BAF and RR are "<<sd_baf<<" "<<sd_rr<<endl;
//  for(int j=0;j<states;++j){
//    sd_baf_vec[j] = sqrt(sd_baf_vec[j]/denom_vec[j]);
//    sd_rr_vec[j] = sqrt(sd_rr_vec[j]/denom_vec[j]);
//    state_object_vector[j].sd_baf = sd_baf_vec[j];
//    state_object_vector[j].sd_rr = sd_rr_vec[j];
//    if (state_object_vector[j].alpha_block>2){
//      //state_object_vector[j].sd_baf = baf_sd; 
//      //state_object_vector[j].sd_rr = rr_sd;
//      //state_object_vector[j].sd_baf = sd_baf;
//      //state_object_vector[j].sd_rr = sd_rr;
//    }
//    //cerr<<"MLE SDs for state "<<j<<" BAF and RR are "<<sd_baf_vec[j]<<" "<<sd_rr_vec[j]<<endl;
//  }
//  //exit(0);
//  //init_emission_cache_mat();
//
//}
//
float HMM_CNV::get_lrr_offset(){
  //return 0;
  float lrr_offset = 0;
  float denom = 0;
  float epsilon = .01;
  for(int i=0;i<total_markers;++i){
    float & lrr = lrr_vec[current_sample*total_markers+i];
    float & baf = baf_vec[current_sample*total_markers+i];
    //float rr = pow(2,lrr);
    //cerr<<"LRR is "<<lrr<<" and mean "<<state_mu_lrr[normal_state]<<endl;
    if (abs(baf-.5)<epsilon){
      float & weight = emission_cache_mat[current_sample*hmm_matrix_size+i*states+normal_state]; 
      lrr_offset+=weight * (lrr-state_mu_lrr[normal_state]);
      denom+=weight; 
    }
  }
  lrr_offset = denom>0?lrr_offset/denom:0;
  cerr<<"Normal state "<<normal_state<<" had mu "<<state_mu_lrr[normal_state]<<" and new LRR offset: "<<lrr_offset<<endl;
  return lrr_offset;
}

//void HMM_CNV::get_alpha_mle(){
//  float denom = 0;
//  map<float,float> backup;
//  for(map<float,float>::iterator it = alpha_weight_map.begin();it!=alpha_weight_map.end();it++){
//    backup[it->first] = it->second;
//    alpha_weight_map[it->first] = 0;
//  }
//  float tally[states];
//  memset(tally,0,sizeof(tally));
//  double SMALL=-1e40;
//  ofstream ofs_alpha("alpha.debug");
//  for(int i=0;i<total_markers;++i){
//    // find the max state
//    int maxindex = -1;
//    double maxprob = SMALL;
//    for(int j=0;j<states;++j){
//      state_t s = state_object_vector[j];
//      //if (emission_cache_mat[i*states+j] > maxprob){
//      if (prob_mat_log[i*states+j] > maxprob){
//        maxindex = j;
//        //maxprob = emission_cache_mat[i*states+j];
//        maxprob = prob_mat_log[i*states+j];
//      }
//    }
//    int max_alpha_block = state_object_vector[maxindex].alpha_block;
//    //state_t s1 = state_object_vector[backtrace_vec[i]];
//    //if (s1.alpha_block>2) cerr<<"Observation:"<<i<<"\t"<<s1.normal_bac<<"\t"<<s1.tumor_cn<<"\t"<<s1.tumor_bac<<"\t"<<s1.alpha<<"\t"<<s1.alpha_block<<endl;
//    //cerr<<"Max alpha block: "<<max_alpha_block<<endl;
//    alpha_block_vec[i] = max_alpha_block;
//    for(int j=0;j<states;++j){
//      state_t s = state_object_vector[j];
//      float alpha = s.alpha;
//      //if (alpha==0.3f){
//        //cerr<<"Adding\n";
//        float weight = (alpha_block_vec[i]>2 && alpha_block_vec[i]==s.alpha_block)?
//        //(emission_cache_mat[i*states+j]):0;
//        exp(prob_mat_log[i*states+j]):0;
//        alpha_weight_map[alpha]+=weight;
//        //if (alpha_block_vec[i]>2 && weight>0){
//          //ofs_alpha<<"Observation:"<<i<<" "<<s.tumor_cn<<" "<<pow(2,lrr_vec[i])<<" "<<s.mu_rr<<" "<<baf_vec[i]<<" "<<s.mu_baf<<" "<<s.alpha<<" "<<weight<<endl;
//        //}
//        tally[j]+=weight;
//      //}else{
//        //cerr<<"alpha: "<<alpha<<endl;
//      //}
//    }
//  }
//  ofs_alpha.close();
//  //exit(0);
//  //cerr<<"Alpha tally:\n";
//  for(int j=0;j<states;++j){
//     //cerr<<j<<"\t"<<tally[j]<<endl;
//  }
//  denom = 0;
////  float beta_a = 1,beta_b=2;
//  for(map<float,float>::iterator it = alpha_weight_map.begin();it!=alpha_weight_map.end();it++){
////    double alpha_prob = gsl_ran_beta_pdf(it->first,beta_a,beta_b);
////    alpha_weight_map[it->first] *= alpha_prob;
//    denom+=alpha_weight_map[it->first];
//  }
//  cerr<<"New alpha prob:";
//  if (denom){
//    cerr<<"updating new\n";
//    mu_alpha = 0;
//    mu_alpha_prob = 0;
//    for(map<float,float>::iterator it = alpha_weight_map.begin();it!=alpha_weight_map.end();it++){
//      alpha_weight_map[it->first]/= denom;
//      //if (alpha_weight_map[it->first] <.001)  alpha_weight_map[it->first] = .001;
//      if (alpha_weight_map[it->first] > mu_alpha_prob) {
//        //mu_alpha = it->first;
//        //mu_alpha_prob = alpha_weight_map[it->first];
//      }
//      mu_alpha+=alpha_weight_map[it->first] * it->first;
//      mu_alpha_prob+= alpha_weight_map[it->first];
//      cerr<<" "<<it->first<<":"<<alpha_weight_map[it->first];
//    }
//    cerr<<endl;
//  }else{
//    cerr<<"rolling back\n";
//    for(map<float,float>::iterator it = alpha_weight_map.begin();it!=alpha_weight_map.end();it++){
//      alpha_weight_map[it->first] = backup[it->first] ;
//      cerr<<" "<<it->first<<":"<<alpha_weight_map[it->first];
//    }
//    cerr<<endl;
//  }
//  //mu_alpha_prob/=alpha_weight_map.size();
//  //cerr<<"size: "<<alpha_weight_map.size()<<" Mean alpha and prob: "<<mu_alpha<<" "<<mu_alpha_prob<<endl;
//  //init_emission_cache_mat();
//}

double unif(){
  return rand()/(RAND_MAX+1.);
}

void HMM_CNV::init_gpu(){
  #ifdef USE_GPU
  cl_int err;
  vector<cl::Platform> platforms;
  err = cl::Platform::get(&platforms);
  // Iterate over platforms
  //ofs<< "Available platforms: " << platforms.size() << endl;
  for(uint i=0;i<platforms.size();++i){
    cerr<<"Platform name "<<platforms[i].getInfo<CL_PLATFORM_NAME>().c_str()<<endl;
    vector<cl::Device> devices;
    platforms[i].getDevices(CL_DEVICE_TYPE_ALL, &devices);
    cerr<<"Available devices: " << devices.size() << endl;
    for  (uint j = 0; j <devices.size(); ++j) {
      string deviceName = devices[j].getInfo<CL_DEVICE_NAME>();
      cl_device_type dtype = devices[j].getInfo<CL_DEVICE_TYPE>();
      switch(dtype){
        case CL_DEVICE_TYPE_GPU:
          cerr<<"GPU "<<j<<"\n";
          break;
        case CL_DEVICE_TYPE_CPU:
          cerr<<"CPU "<<j<<"\n";
          break;
      }
    }
  }
   // create a context from specified platform
  cl_context_properties cps[3] = {CL_CONTEXT_PLATFORM,(cl_context_properties)(platforms[platform_id])(),0};
  cl::Context context;
  cl_device_type  dtype = platform_type.compare("cpu")?CL_DEVICE_TYPE_GPU:CL_DEVICE_TYPE_CPU;
  //context = platform_id==0?cl::Context(CL_DEVICE_TYPE_CPU,cps):cl::Context(CL_DEVICE_TYPE_GPU,cps);
  context = cl::Context(dtype,cps);
  cerr<<"Selected platform ID: "<<platform_id<<endl;
  // extract the devices from this platform
  vector<cl::Device> devices = context.getInfo<CL_CONTEXT_DEVICES>();
  int device_id = 0;
  // make a command queue from the first device
  opencl_info.command_queue = cl::CommandQueue(context,devices[device_id]);
  // Create a program from the kernel source
  ifstream ifs(kernel_path.data());
  clSafe(ifs.is_open()?CL_SUCCESS:-1,"kernel_path not found");
  string source_str(istreambuf_iterator<char>(ifs),(istreambuf_iterator<char>()));
  cl::Program::Sources source(1,make_pair(source_str.c_str(),source_str.length()+1));
  cl::Program program(context,source);
  err = program.build(devices);
  cerr<<"Building kernel\n";
  if (err!=CL_SUCCESS) {
      cerr<<"Build failed:\n";
      string buffer;
      program.getBuildInfo(devices[device_id],CL_PROGRAM_BUILD_LOG,&buffer);
      cerr<<buffer<<endl;
  }else{
    cerr<<"Built kernel\n";
  }
 // Create the OpenCL kernels
  opencl_info.kernel_rescale_transition = cl::Kernel(program, "rescale_transition", &err);
  clSafe(err,"created kernel");
  opencl_info.kernel_forward_backward = cl::Kernel(program, "forward_backward", &err);
  clSafe(err,"created kernel");
  //opencl_info.kernel_viterbi = cl::Kernel(program, "viterbi", &err);
  //clSafe(err,"created kernel");
  opencl_info.kernel_likelihood = cl::Kernel(program, "likelihood", &err);
  clSafe(err,"created kernel");
  opencl_info.kernel_posterior = cl::Kernel(program, "posterior", &err);
  clSafe(err,"created kernel");
  opencl_info.kernel_baumwelch_normalize = cl::Kernel(program, "baumwelch_normalize", &err);
  clSafe(err,"created kernel");
  opencl_info.kernel_baumwelch_transition = cl::Kernel(program, "baumwelch_transition", &err);
  clSafe(err,"created kernel");
  opencl_info.kernel_zero_transition = cl::Kernel(program, "zero_transition", &err);
  clSafe(err,"created kernel");
  opencl_info.kernel_normalize_transition = cl::Kernel(program, "normalize_transition", &err);
  clSafe(err,"created kernel");
// Create the memory buffers
  opencl_info.cn_vec_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (states)*sizeof(int), NULL, &err);
  clSafe(err, "created buffer for cn vec");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.cn_vec_buffer, CL_TRUE, 0,  sizeof(int)*(states), state_tumorcn_vec , NULL, NULL );
  clSafe(err, "wrote buffer for cn vec");
  opencl_info.penalty_mat_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (states*states)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for penalty matrix");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.penalty_mat_buffer, CL_TRUE, 0,  sizeof(float)*(states*states), penalty_mat , NULL, NULL );
  clSafe(err, "wrote buffer for penalty matrix");
  opencl_info.transition_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states*states)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for transition matrix");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.transition_mat_buffer, CL_TRUE, 0,  sizeof(float)*(states*states), trans_mat , NULL, NULL );
  clSafe(err, "wrote buffer for transition matrix");
  opencl_info.transition_new_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states*states)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for new transition matrix");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.transition_new_mat_buffer, CL_TRUE, 0,  sizeof(float)*(states*states), trans_mat , NULL, NULL );
  clSafe(err, "wrote buffer for new transition matrix");
  opencl_info.scaled_forward_transition_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states*states)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for transition matrix");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.scaled_forward_transition_mat_buffer, CL_TRUE, 0,  sizeof(float)*(states*states), trans_mat , NULL, NULL );
  clSafe(err, "wrote buffer for transition matrix");
  opencl_info.scaled_backward_transition_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states*states)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for transition matrix");
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.scaled_backward_transition_mat_buffer, CL_TRUE, 0,  sizeof(float)*(states*states), trans_mat , NULL, NULL );
  clSafe(err, "wrote buffer for transition matrix");
  //opencl_info.emission_mat_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (states*symbols)*sizeof(float), NULL, &err);
  //clSafe(err, "created buffer for emission matrix");
  //opencl_info.emission_mat_new_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states*symbols)*sizeof(float), NULL, &err);
  //clSafe(err, "created buffer for new emission matrix");
  //err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.emission_mat_buffer, CL_TRUE, 0,  sizeof(float)*(states*symbols), emission_mat , NULL, NULL );
  //clSafe(err, "wrote buffer for emission matrix");
  opencl_info.emission_cache_mat_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (states*total_markers)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for emission cache matrix");
  opencl_info.forward_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_markers*states)*sizeof(double), NULL, &err);
  clSafe(err, "created buffer for forward matrix");
  //err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_mat_buffer, CL_TRUE, 0,  sizeof(double)*(states*total_markers), forward_mat_log_t , NULL, NULL );
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.forward_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), forward_mat_log , NULL, NULL );
  clSafe(err, "wrote buffer for forward matrix");
  opencl_info.backward_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_markers*states)*sizeof(double), NULL, &err);
  clSafe(err, "created buffer for backward matrix");
  //err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_mat_buffer, CL_TRUE, 0,  sizeof(double)*(states*total_markers), backward_mat_log_t , NULL, NULL );
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.backward_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), backward_mat_log , NULL, NULL );
  clSafe(err, "wrote buffer for backward matrix");
  opencl_info.prob_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_markers*states)*sizeof(double), NULL, &err);
  clSafe(err, "created buffer for prob matrix");
//err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.prob_mat_buffer, CL_TRUE, 0,  sizeof(double)*(states*total_markers), prob_mat_log_t , NULL, NULL );
  err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.prob_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), prob_mat_log , NULL, NULL );
  clSafe(err, "wrote buffer for prob matrix");
  //opencl_info.greedy_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_markers*states)*sizeof(double), NULL, &err);
  //clSafe(err, "created buffer for greedy matrix");
  //err = opencl_info.command_queue.enqueueWriteBuffer(opencl_info.greedy_mat_buffer, CL_TRUE, 0,  sizeof(double)*(total_markers*states), greedy_mat_log , NULL, NULL );
  //clSafe(err, "wrote buffer for greedy matrix");
  //opencl_info.bestpath_mat_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_markers*states)*sizeof(int), NULL, &err);
  //clSafe(err, "created buffer for bestpath matrix");
  //opencl_info.backtrace_vec_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (total_markers)*sizeof(int), NULL, &err);
  //clSafe(err, "created buffer for backtrace vec");
  opencl_info.obs_index_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (1*1)*sizeof(int), NULL, &err);
  clSafe(err, "created buffer for observation index");
  opencl_info.counter_index_buffer = cl::Buffer(context, CL_MEM_READ_ONLY, (states*1)*sizeof(int), NULL, &err);
  clSafe(err, "created buffer for counter index");
  opencl_info.likelihood_log_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (1*1)*sizeof(float), NULL, &err);
  clSafe(err, "created buffer for likelihood log");
  opencl_info.baumwelch_norm_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, (states)*sizeof(double), NULL, &err);
  clSafe(err, "created buffer for Baum Welch norm");
  opencl_info.forward_gap_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, sizeof(float), NULL, &err);
  clSafe(err, "created buffer for gap forward");
  opencl_info.backward_gap_buffer = cl::Buffer(context, CL_MEM_READ_WRITE, sizeof(float), NULL, &err);
  clSafe(err, "created buffer for gap backward");

  // CREATE LOCAL BUFFERS

  int arg;
  arg = 0;
  err = opencl_info.kernel_rescale_transition.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.cn_vec_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.forward_gap_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.backward_gap_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.penalty_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.scaled_forward_transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, opencl_info.scaled_backward_transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, cl::__local(sizeof(float)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_rescale_transition.setArg(arg++, cl::__local(sizeof(float)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");

  arg = 0;
  err = opencl_info.kernel_forward_backward.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, total_markers);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, opencl_info.counter_index_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, opencl_info.scaled_forward_transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, opencl_info.scaled_backward_transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, opencl_info.emission_cache_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, opencl_info.forward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, opencl_info.backward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_forward_backward.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");

  //arg = 0;
//  err = opencl_info.kernel_viterbi.setArg(arg++, states);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, total_markers);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, opencl_info.counter_index_buffer);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, opencl_info.scaled_forward_transition_mat_buffer);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, opencl_info.emission_cache_mat_buffer);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, opencl_info.greedy_mat_buffer);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, opencl_info.bestpath_mat_buffer);
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
//  clSafe(err,"clSetKernelArg");
//  err = opencl_info.kernel_viterbi.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
//  clSafe(err,"clSetKernelArg");

  arg = 0;
  err = opencl_info.kernel_likelihood.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_likelihood.setArg(arg++, opencl_info.forward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_likelihood.setArg(arg++, opencl_info.backward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_likelihood.setArg(arg++, opencl_info.likelihood_log_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_likelihood.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_likelihood.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_likelihood.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");

  arg = 0;
  cerr<<"Kernel arguments posterior\n";
  err = opencl_info.kernel_posterior.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_posterior.setArg(arg++, opencl_info.forward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_posterior.setArg(arg++, opencl_info.backward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_posterior.setArg(arg++, opencl_info.likelihood_log_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_posterior.setArg(arg++, opencl_info.prob_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_posterior.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_posterior.setArg(arg++, cl::__local(sizeof(double)*HMM_STATE_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  arg = 0;
  cerr<<"Kernel arguments baumwelch normalize\n";
  err = opencl_info.kernel_baumwelch_normalize.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_normalize.setArg(arg++, total_markers);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_normalize.setArg(arg++, opencl_info.prob_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_normalize.setArg(arg++, opencl_info.baumwelch_norm_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_normalize.setArg(arg++, cl::__local(sizeof(double)*HMM_OBS_BLOCK_WIDTH));
  clSafe(err,"clSetKernelArg");
  arg = 0;
  err = opencl_info.kernel_zero_transition.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_zero_transition.setArg(arg++, opencl_info.transition_new_mat_buffer);
  clSafe(err,"clSetKernelArg");
  arg = 0;
  cerr<<"Kernel arguments transition\n";
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.counter_index_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.forward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.backward_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.emission_cache_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.likelihood_log_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.transition_new_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_baumwelch_transition.setArg(arg++, opencl_info.scaled_forward_transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  arg = 0;
  err = opencl_info.kernel_normalize_transition.setArg(arg++, states);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_normalize_transition.setArg(arg++, opencl_info.transition_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_normalize_transition.setArg(arg++, opencl_info.transition_new_mat_buffer);
  clSafe(err,"clSetKernelArg");
  err = opencl_info.kernel_normalize_transition.setArg(arg++, opencl_info.baumwelch_norm_buffer);
  clSafe(err,"clSetKernelArg");
  #endif
}

int main_hmm(int argc, char * argv[]){
  return 0;
}
